# This file takes the individual level output and formats it
# The individual level results build on data that cannot be made available for replication. 

load("ind_effect_all.rdata")
load("year_effect_all.rdata")

# This file esimates the effect conditional on cohabitation using only municipalities with recurring data

#  install.packages("xtable")

library(xtable)

# create table with findings

main_tab <- matrix(NA, ncol = 10, nrow = 11)

# first input yearly quantities
for(i in 1:4){
  main_tab[1 , c(i, i + 5)] <- round(year_effect[c(4 * (i - 1) + 1, 4 * (i - 1) + 3), 1], 4) * 100
  main_tab[2 , c(i, i + 5)] <- round(year_effect[c(4 * (i - 1) + 1, 4 * (i - 1) + 3), 2], 4) * 100
  main_tab[4 , c(i, i + 5)] <- round(year_effect[c(4 * (i - 1) + 2, 4 * (i - 1) + 4), 1], 4) * 100
  main_tab[5 , c(i, i + 5)] <- round(year_effect[c(4 * (i - 1) + 2, 4 * (i - 1) + 4), 2], 4) * 100
  main_tab[7 , c(i, i + 5)] <- round(year_effect[c(4 * (i - 1) + 2, 4 * (i - 1) + 4), 6], 4) * 100
}

# then input pooled quantities

main_tab[8:9, c(5, 10)] <- c(sum(main_tab[8:9, 1:4]))

main_tab[7, c(5, 10)] <- round(main_effect[, 5], 4) * 100

main_tab[1, c(5, 10)] <- round(main_effect[, 1], 4) * 100
main_tab[2, c(5, 10)] <- round(main_effect[, 2], 4) * 100

main_tab[4, c(5, 10)] <- round(main_effect[, 3], 4) * 100
main_tab[5, c(5, 10)] <- round(main_effect[, 4], 4) * 100

main_tab[c(3, 6), ] <- rbind(paste("[",  round(qnorm(0.025) * main_tab[2, ] + main_tab[1, ], 2) ,
                                   ", ", round(qnorm(0.975) * main_tab[2, ] + main_tab[1, ], 2) ,
                                   "]", sep = ""),
                             paste("[",  round(qnorm(0.025) * main_tab[5, ] + main_tab[4, ], 2),
                                   ", ", round(qnorm(0.975) * main_tab[5, ] + main_tab[4, ], 2),
                                   "]", sep = ""))

main_tab[7, c(5, 10)] <- round(main_effect[, 5], 4) * 100

for(i in 1:4){
  main_tab[8 , c(i, i + 5)] <- paste(year_effect[c(4 * (i - 1) + 1, 4 * (i - 1) + 3), 3])
  main_tab[9 , c(i, i + 5)] <- paste(year_effect[c(4 * (i - 1) + 1, 4 * (i - 1) + 3), 4])
  main_tab[10, c(i, i + 5)] <- paste(round(year_effect[c(4 * (i - 1) + 1, 4 * (i - 1) + 3), 5]) )
  main_tab[11, c(i, i + 5)] <- paste(round(year_effect[c(4 * (i - 1) + 2, 4 * (i - 1) + 4), 5]) )
  
}

main_tab[8, c(5, 10)] <- c(paste(sum(year_effect[seq(1, 13, 4), 3])), paste(sum(year_effect[seq(3, 15, 4), 3])))
main_tab[9, c(5, 10)] <- c(paste(sum(year_effect[seq(1, 13, 4), 4])), paste(sum(year_effect[seq(3, 15, 4), 4])))

rownames(main_tab) <- c("Estimate", "St.err.", "CI",
                        "Robust estimate", "Robust st.err.", "Robust CI",
                        "Mean turnout among parents of ineligible", "N_ineligible",
                        "N_eligible", "Bandwidth", "Bias-correction bandwidth")

colnames(main_tab) <- rep(c("2009", "2013", "2014", "2015", "Pooled Effect"), 2)

xtable(main_tab)

